Проблемой короновируса COVID-19 буквально сейчас заняты многие ученые по всему миру. В стороне не осталось и комьюнити Data scientist-ов, которые пытаются не просто предсказать развитие пандемии, а помочь определить наиболее значимые факторы, влияющие на распространение заразы. Например на небезызвестном Kaggle.
В данной домашней работе Вы будете использовать ежедневно обновляемые данные из репозитория, а конкретно time-series-19-covid-combined.csv.
Вам предстоит изучить имеющиеся данные, попробовать найти некоторые закономерности и ответить на вопросы. Ниже представлены задания как исследовательского характера, где приветствуется инициативность, так и чисто тренировочного харатера - где необходимо использовать изученные на семинаре инструменты для работы с временными рядами. Не забывайте сопровождать ваш код развернутыми комментариями и выводами. Чем их больше и чем они качественнее, тем лучше и выше вероятность получения максимальной оценки за домашнее задание. Не болейте!
!wget -O covid.csv https://raw.githubusercontent.com/datasets/covid-19/master/data/time-series-19-covid-combined.csv --no-check-certificate -q
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
pd.options.display.float_format = '{:.2f}'.format
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (15, 15)
df = pd.read_csv('covid.csv', parse_dates=['Date'])
df.head()
Границы исследования
df['Date'].agg(['min', 'max'])
Количество заболевших, выздоровивших и погибших
df[df['Date'] == df['Date'].max()][['Confirmed', 'Recovered', 'Deaths']].sum()
Карта
Данные логарифмированы, иначе будет одна яркая точка -- США. На карте доступна анимация, а так же масштабирвоание.
import folium
import folium.plugins as plugins
data = []
__df = df[(~df['Confirmed'].isna()) & (df['Province/State'] != 'Diamond Princess')].copy()
__df['Confirmed'] = np.log(1 + __df['Confirmed'])
__df['Confirmed'] = __df['Confirmed'] / __df['Confirmed'].max()
__df = __df[['Date', 'Lat', 'Long', 'Confirmed']].groupby('Date').agg(lambda x: list(x))
for date, row in __df.iterrows():
row = list(zip(*row))
row = [[x[0], x[1], np.clip(x[2], 0.00001, 1)] for x in row if x[2] > 0]
data.append(row)
m = folium.Map([0., 0.], zoom_start=2)
hm = plugins.HeatMapWithTime(data,
index=list(__df.index.strftime('%Y-%m-%d').values),
auto_play=True,
min_speed=1,
)
hm.add_to(m)
m
Примечание. В анализе временных рядов иногда необходимо изменить частоту дискретизации (Resampling). Такая операция может понадобится, когда имеющихся данных недостаточно или их, наоборот, слишком много. Может не устраивать имеющаяся частота или просто хочется посмотреть на данные с более общего ракурса.
Различают два вида изменения частоты: повышение (Upsampling) и понижение (Downsampling). При повышении временной ряд пересчитывается с низкой частоы на более высокую частоту (например от годовой до месячной частоты). В таком случае процесс включает в себя заполнение или интерполяцию появившихся пропусков в данных. При понижении временной ряд передискретизируется с высокой частоты на низкую (наример с еженедельной на месячную частоту). Это включает в себя агрегацию существующих данных.
data = df[df['Confirmed'] > 0].groupby('Date')['Country/Region'].nunique()
data.name = 'Approved countries'
data.index = data.index.strftime('%Y-%m-%d')
data.plot.bar(figsize=(15, 7), title='Количество зараженных стран по дням');
data = df[df['Confirmed'] > 0].groupby(pd.Grouper(key='Date', freq='W-MON'))['Country/Region'].nunique()
data.name = 'Approved countries'
data.index = data.index.strftime('%Y-%m-%d')
data.plot.bar(figsize=(15, 7), title='Количество зараженных стран по неделям');
data = df.groupby('Date')[['Confirmed', 'Deaths', 'Recovered']].sum()
data.index = data.index.strftime('%Y-%m-%d')
data.plot(figsize=(15, 7), title='Количество заболевших | погибших | выздоровивших');
__df = df[(df['Date'] == df['Date'].max()) & (df['Confirmed'] > 0)].copy()
__df = __df.groupby('Country/Region')[['Deaths', 'Confirmed']].sum()
__df['cfr'] = __df['Deaths'] / __df['Confirmed']
__df.sort_values('cfr', inplace=True)
__df = __df[__df['Confirmed'] > 100]
# __df = __df[__df['cfr'] > __df['cfr'].mean()]
__df['cfr'].plot.barh(figsize=(7, 30), title='Летальность в разных странах');
Подчистим данные
df_c = df.groupby(['Date', 'Country/Region'])[['Confirmed', 'Recovered', 'Deaths']].sum().reset_index()
df_c
region_start = df_c[df_c['Confirmed'] > 0].groupby('Country/Region')['Date'].min()
region_start.name = 'd_start'
__df = pd.merge(df_c[df_c['Confirmed'] > 0], region_start, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start']).dt.days
__df['Active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Confirmed', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Confirmed']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Confirmed_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row['Confirmed'], label=f'{i+1:>2}.{name:<15}- total {row["Confirmed_max"]:>8,d}' if i < 20 else None)[0].get_color()
if i < 10:
plt.scatter(row['n_day'][-1], row['Confirmed'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Confirmed'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Число зараженных')
plt.title('Распространение болезни в разных странах с первого подтвержденного случая');
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Active', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Active']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Active_max', ascending=False)
).iterrows()):
if i < 10:
color = plt.plot(row['n_day'], row['Active'], label=f'{i+1:>2}.{name:<15}- total {row["Active_max"]:>8,d}' if i < 20 else None)[0].get_color()
plt.scatter(row['n_day'][-1], row['Active'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Active'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Число активных больных')
plt.title('Число активно больных в разных странах с первого подтвержденного случая');
(
__df
.groupby('Country/Region')['n_day']
.max()
.sort_values(ascending=False)
.plot
.barh(figsize=(7, 30), title='Сколько дней прошло с первого заражения')
);
region_start_100 = df_c[df_c['Confirmed'] >= 100].groupby('Country/Region')['Date'].min()
region_start_100.name = 'd_start100'
__df = pd.merge(df_c[df_c['Confirmed'] >= 100], region_start_100, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start100']).dt.days
__df['Active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Confirmed', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Confirmed']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Confirmed_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row['Confirmed'], label=f'{i+1:>2}.{name:<15}- total {row["Confirmed_max"]:>8,d}' if i < 20 else None)[0].get_color()
if i < 10:
plt.scatter(row['n_day'][-1], row['Confirmed'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Confirmed'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel('Число зараженных')
plt.title('Распространение болезни в разных странах с 100 подтвержденного случая');
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Active', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Active']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Active_max', ascending=False)
).iterrows()):
if i < 10:
color = plt.plot(row['n_day'], row['Active'], label=f'{i+1:>2}.{name:<15}- total {row["Active_max"]:>8,d}' if i < 20 else None)[0].get_color()
plt.scatter(row['n_day'][-1], row['Active'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Active'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Число активных больных')
plt.title('Число активно больных в разных странах с первого подтвержденного случая');
(
__df
.groupby('Country/Region')['n_day']
.max()
.sort_values(ascending=False)
.plot
.barh(figsize=(7, 30), title='Сколько дней прошло с 100 заражения')
);
Особой разницы между графиками не заметно, но от 100 зараженного выглядит более гладким
Если смотреть на Китай и Германию, то на переход от начала заражения до пика проходит меньше времени, чем от пика до спада.
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['Confirmed', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['Confirmed']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('Confirmed_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row['Confirmed'], label=f'{i+1:>2}.{name:<15}- total {row["Confirmed_max"]:>8,d}' if i < 20 else None)[0].get_color()
if 1 <= i < 15:
plt.scatter(row['n_day'][-1], row['Confirmed'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['Confirmed'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel('Число зараженных')
plt.ylim([-4000, 1.1 * __df.groupby('Country/Region')['Confirmed'].max().sort_values()[-2]])
plt.title('Распространение болезни в разных странах с 100 подтвержденного случая');
plt.tight_layout()
Количество подтвержденных случаев в России близко к аналогичным в Бразилии, Канаде, Китаю и Ирану.
Если сравнивать темпы роста, то Россия похожа на UK, Францию, Италию, Германию.
По предыдущим графикам, у России еще все впереди.
for col in ['Confirmed', 'Recovered', 'Deaths', 'Active']:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[[col, 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')[col]
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values(f'{col}_max', ascending=False)
.head(10)
).iterrows()):
color = plt.plot(row['n_day'], row[col], label=f'{i+1:>2}.{name:<15}- total {row[f"{col}_max"]:>8,d}' if i < 20 else None)[0].get_color()
if 0 <= i < 10:
plt.scatter(row['n_day'][-1], row[col][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row[col][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel(f'Число {col}')
# plt.ylim([-4000, 1.1 * __df.groupby('Country/Region')[col].max().sort_values()[-2]])
plt.title(f'Распространение болезни ({col}) в разных странах с 100 подтвержденного случая');
plt.tight_layout()
plt.show()
Все показатели, кроме выздоровивших похожи. Лидирует США. Во некоторых странах ситуация намного хуже чем в Китае.
По количеству выздоровивших неожиданно лидирует Германия.
import seaborn as sns
def plot_corr(D, size):
corr = D.corr()
# corr = np.abs(corr)
f, ax = plt.subplots(figsize=(size, size))
cmap = plt.cm.Oranges
sns.heatmap(corr, cmap=cmap,
xticklabels=corr.columns,
yticklabels=corr.columns)
from IPython.display import display
for col in ['Confirmed', 'Recovered', 'Deaths', 'Active']:
print(f'Корреляция по признаку {col}')
display(
(__df
.groupby('Country/Region')[[col, 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')[col]
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values(f'{col}_max', ascending=False)
.head(10)
)
.T
.apply(lambda x: pd.Series(x.loc[col]))
.corr()
.style
.background_gradient(cmap='coolwarm')
.set_precision(2)
)
По всем показателям, кроме Recovered Китай не похож на остальные страны. В топе стран высокая кореляция по наблюдаемым величинам. По Recovered все страны похожи
1.7 Для первых 5 стран из топ-10 и России постойте сравнительные графики. Изучите как изменялись значения отношений погибших/выздоровивших, погибших/заболевших, заболевших/общее число жителей, ваш вариант...
!wget -O population.csv https://raw.githubusercontent.com/datasets/population/master/data/population.csv --no-check-certificate -q
top_c = list(set(df_c.groupby('Country/Region')['Confirmed'].max().sort_values(ascending=False).head(5).index) | {'Russia'})
top_c
population = pd.read_csv('population.csv', sep=',')
population['Country Name'].replace({'United States': 'US', 'Russian Federation': 'Russia', 'Iran, Islamic Rep.': 'Iran'}, inplace=True)
population = population[population['Country Name'].isin(top_c)]
population = population.loc[population.groupby('Country Name')['Year'].idxmax().values].set_index('Country Name')['Value']
population.name = 'population'
population.index.name = 'Country/Region'
set(population.index) - set(top_c)
__df = df_c[df_c['Country/Region'].isin(top_c)]
# region_start_100 = df_c[df_c['Confirmed'] >= 100].groupby('Country/Region')['Date'].min()
# region_start_100.name = 'd_start100'
__df = __df.join(region_start_100, on='Country/Region').join(population, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start100']).dt.days
__df = __df[__df['n_day'] > 0]
__df['Active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df['Deaths/Recovered'] = __df['Deaths'] / (__df['Recovered'] + 1) # prevent divide by 0
__df['Deaths/Confirmed'] = __df['Deaths'] / __df['Confirmed']
__df['Confirmed/population'] = __df['Confirmed'] / __df['population']
__df.head()
for col in ['Deaths/Recovered', 'Deaths/Confirmed', 'Confirmed/population', 'Active']:
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[[col, 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')[col]
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values(f'{col}_max', ascending=False)
).iterrows()):
color = plt.plot(row['n_day'], row[col], label=f'{i+1:>2}.{name:<10}')[0].get_color()
plt.scatter(row['n_day'][-1], row[col][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row[col][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 100 зараженного')
plt.ylabel(f'Число {col}')
plt.title(f'Распространение болезни ({col}) в разных странах с 100 подтвержденного случая');
plt.tight_layout()
plt.show()
Deaths/Recovered: С начала заболевания много людей погибает, а количество выздоровивших растет с некоторой задержкой, поэтому есть всплеск данного показателя в начале болезни в странах.
Deaths/Confirmed: В России подозрительно низкий показатель смертности. Если построить аналогичный график для 25 стран, то у России данный показатель также низок, на уровне Саудовской Аравии (в отличии от России, у нее он падал со временем). По всей видимости, в России по другому считают стасистику погибших (умерших от вируса, а не умерших с вирусом). Также, возможно, занижают показатели для более низкого значения смертности
Confirmed/population: Видно, что в США хоть и больше всего зараженных, но население страны так же большое. По данному показателю лидирует Испания. Также заметно, что заболело в каждой стране менее 1% населения, а в России менее 0,1%.
Active: В Германии заболевание идет на спад, в Италии, Испании и Франции рост прекращается. По России трудно говорить о начале пика заболевших. В США все плохо.
1.8 Выделите временной ряд по одному из целевых признаков. Для выделенного временного ряда:
__df = pd.concat([df_c[['Date', 'Country/Region', 'Confirmed']], df_c.shift(df_c.groupby('Date').size().unique()[0])['Confirmed']], axis=1).dropna()
__df.columns = ['Date', 'Country/Region', 'Confirmed', 'Confirmed_lag']
__df['lag'] = __df['Confirmed'] - __df['Confirmed_lag']
region_start = __df[__df['Confirmed'] > 0].groupby('Country/Region')['Date'].min()
region_start.name = 'd_start'
__df = pd.merge(__df[__df['Confirmed'] > 0], region_start, on='Country/Region')
__df['n_day'] = (__df['Date'] - __df['d_start']).dt.days
__df['day_of_week'] = __df['Date'].dt.dayofweek
__df
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df
.groupby('Country/Region')[['lag', 'n_day']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['lag']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('lag_max', ascending=False)
.head(10)
).iterrows()):
color = plt.plot(row['n_day'], row['lag'], label=f'{i+1:>2}.{name:<15}- total {row["lag_max"]:>8,d}' if i < 20 else None)[0].get_color()
plt.scatter(row['n_day'][-1], row['lag'][-1], color=color, alpha=1)
plt.annotate(f'{name}', (row['n_day'][-1], row['lag'][-1]),
xytext=(3, 7), textcoords='offset points', ha='center', va='bottom', color=color,
bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.2))
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с первого зараженного')
plt.ylabel('Прирост больных')
plt.title('Число прироста новых случаев в разных странах');
Ряд нестационарен, заметны некоторые циклические колебания (по дням недели?).
Теоретически не может быть стационарным, когда-то все выздоровят. В Китае количество новых случаев давно стало нулевым
1.9 Для любого интересующего Вас временного ряда постройте графики ACF и PACF. Сделайте выводы.
from statsmodels.tsa.stattools import acf, pacf
pd.Series(acf(__df[__df['Country/Region'] == 'Russia']['lag'])).plot.bar(figsize=(14,7), title='ACF');
pd.Series(pacf(__df[__df['Country/Region'] == 'Russia']['lag'])).plot.bar(figsize=(14,7), title='PACF');
оптимальный параметры для модели ARMA (1, 2)
Примеры признаков - количество дней со начала карантина, конинент, плотность населения, количество употребляемого алкоголя, летучих мышей, алкоголя, крокодилов на душу населения, количество туристов в год, уровень безработицы, температура, количество игроков в Plague Inc., средний возраст или продолжительность жизни и так далее. Вы можете использовать любой найденный датасет. Главное, чтобы у Вас получилось скомбинировать данные.
Примеры датасеты:
population = pd.read_csv('population.csv', sep=',')
population['Country Name'].replace(
{
'United States': 'US',
'Russian Federation': 'Russia',
'Iran, Islamic Rep.': 'Iran',
'Czech Republic': 'Czechia',
'Egypt, Arab Rep.': 'Egypt',
'Korea, Rep.': 'Korea, South',
'Slovak Republic': 'Slovakia',
}, inplace=True)
population = population.loc[population.groupby('Country Name')['Year'].idxmax().values].set_index('Country Name')['Value']
population.name = 'population'
population.index.name = 'Country/Region'
__df = df_c.copy()
__df['cfr'] = __df['Deaths'] / __df['Confirmed']
__df = pd.concat([__df, __df.shift(df_c.groupby('Date').size().unique()[0])['cfr']], axis=1).dropna()
cols = list(__df.columns)
cols[-1] = cols[-1] + '_lag'
__df.columns = cols
__df['lag'] = __df['cfr'] - __df['cfr_lag']
region_start1 = __df[__df['Confirmed'] > 0].groupby('Country/Region')['Date'].min()
region_start1.name = 'd_start1'
region_start100 = __df[__df['Confirmed'] > 100].groupby('Country/Region')['Date'].min()
region_start100.name = 'd_start100'
region_start1000 = __df[__df['Confirmed'] > 1000].groupby('Country/Region')['Date'].min()
region_start1000.name = 'd_start1000'
__df = (
__df[__df['Confirmed'] > 0]
.join(region_start1, on='Country/Region', how='inner')
.join(region_start100, on='Country/Region', how='inner')
.join(region_start1000, on='Country/Region', how='inner')
.join(population, on='Country/Region', how='inner')
)
__df['active'] = __df['Confirmed'] - __df['Recovered'] - __df['Deaths']
__df['n_day1'] = (__df['Date'] - __df['d_start1']).dt.days
__df['n_day100'] = (__df['Date'] - __df['d_start100']).dt.days
__df['n_day1000'] = (__df['Date'] - __df['d_start1000']).dt.days
__df['confirmed_rate'] = __df['Confirmed'] / __df['population']
__df['recovered_rate'] = __df['Recovered'] / __df['population']
__df['deaths_rate'] = __df['Deaths'] / __df['population']
__df['active_rate'] = __df['active'] / __df['population']
__df['day_of_week'] = __df['Date'].dt.dayofweek
__df = __df[__df['n_day100'] >= 0]
__df.drop(columns=['cfr_lag', 'd_start1', 'd_start100', 'd_start1000', 'population'], inplace=True)
__df.reset_index(drop=True, inplace=True)
__df.to_csv('covid_features.csv')
__df.to_pickle('covid_features.pkl')
__df
plt.figure(figsize=(20, 20), facecolor='w')
plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05, hspace=.01)
plot_num = 1
for i, (country, active) in enumerate(__df.groupby('Country/Region')['active'].max().sort_values(ascending=False).head(10).iteritems()):
___df = __df[__df['Country/Region'] == country]
for param in ['n_day1', 'n_day100', 'n_day1000', 'confirmed_rate', 'recovered_rate', 'deaths_rate', 'active_rate', 'day_of_week']:
plt.subplot(10, 8, plot_num)
if i == 0:
plt.title(param, size=18)
if param == 'n_day1':
plt.ylabel(country, fontsize=18)
if ___df.dtypes[param] == np.int64:
___df.groupby(param)['lag'].mean().plot.bar()
else:
___df.groupby(param)['lag'].mean().plot()
plt.xticks(())
plt.yticks(())
plot_num += 1
На графиках представлена зависимость прироста смертности от параметра дня в разных странах
выводы были после каждого пункта
Теперь вам предстоит построить несколько моделей и спрогнозировать временной ряд для такого показателя как летальности от вируса. Летальность можно рассчитывать по формуле:
$$ CFR = \frac{Deaths}{Confirmed} $$Для каждого типа модели сделайте несколько прогнозов для не менее 3 разных стран/регионов/другое (на ваш обоснованный выбор). Главное, чтобы каждый из типов моделей прознозировал одинаковый набор данных, чтобы в конце можно было сделать выводы о качестве работы той или иной модели.
Метрика качества RMSE. Не забудьте разбить данные временной ряд на данные для обучения и тестовые.
Не забывайте сопровождать ваш код комментариями, графиками и выводами.
import datetime
from sklearn.metrics import mean_squared_error
def rmse(y_true, y_pred):
return np.sqrt(mean_squared_error(y_true, y_pred))
plt.figure(figsize=(20,10), facecolor='w')
for i, (name, row) in enumerate(( __df[__df['n_day1000'] > 0]
.groupby('Country/Region')[['lag', 'n_day1000']]
.agg(lambda x: list(x))
.join(
__df.groupby('Country/Region')['n_day1000']
.max()
.astype(np.int32),
rsuffix='_max')
.sort_values('n_day1000_max', ascending=False)
.head(10)
).iterrows()):
color = plt.plot(row['n_day1000'], row['lag'], label=f'{i+1:>2}.{name:<15}')[0].get_color()
plt.scatter(row['n_day1000'][-1], row['lag'][-1], color=color, alpha=1)
plt.legend(prop={'family': 'monospace'})
plt.xlabel('Дней с 1000 зараженного')
plt.ylabel('Прирост летальности')
plt.title('Число прироста летальности в разных странах');
split_date = __df['Date'].max() - datetime.timedelta(10)
print(split_date)
train = __df[__df['Date'] < split_date].copy()
test = __df[__df['Date'] >= split_date].copy()
all_countries = list(train['Country/Region'].unique())
top_countries = ['US', 'Russia', 'Spain']
def measure_quality(df, name):
per_country = (
df
.groupby('Country/Region')
.agg(lambda x: list(x))
.apply(lambda x: rmse(x['cfr'], x['cfr_r']) if np.isnan(x['cfr_r']).sum() == 0 else np.nan, axis=1)
)
per_country.name = 'rmse'
world = pd.merge(per_country, population, left_index=True, right_index=True, how='inner')
per_country['World'] = (world['rmse'] * world['population']).sum() / world['population'].sum()
per_country.name = name
return per_country
class Experimens:
def __init__(self, top_cols):
self._df = pd.DataFrame(columns=['World'] + top_cols)
def add(self, data, name):
quality = measure_quality(data, name)
self._df = self._df.append(quality, verify_integrity=True)
return self.show
@property
def show(self):
return self._df.copy()
ex = Experimens(top_countries)
def plot_predict(predict, title=None):
plt.figure(figsize=(14, 7))
for country in top_countries:
_tr = train[train['Country/Region'] == country]
_te = predict[predict['Country/Region'] == country]
color = plt.plot(_tr['Date'], _tr['cfr'], label=country + ' train')[0].get_color()
plt.plot(_te['Date'], _te['cfr'], label=country + ' test', color=color, marker='+')
plt.plot(_te['Date'], _te['cfr_r'], label=country + ' predict', color=color, marker='o')
plt.title(title)
plt.legend();
predict = train.groupby('Country/Region')['cfr'].mean()
predict = test[['Date', 'Country/Region', 'lag', 'cfr']].join(predict, on='Country/Region', rsuffix='_r')
plot_predict(predict, 'Константное предсказание')
ex.add(predict, 'const')
def flatten(df):
res = []
for index, row in df.iterrows():
for r2 in zip(*row):
new_row = {
df.index.name: index
}
for keys, values in zip(row.keys(), r2):
new_row[keys] = values
res.append(new_row)
return pd.DataFrame(res)
def predict_to_normal(data, data_start):
data_start = data_start[data_start['Date'] == data_start['Date'].max()].set_index('Country/Region')[['cfr']]
df = data[['Date', 'Country/Region', 'lag', 'lag_r']].groupby('Country/Region').agg(lambda x: list(x))
df['cfr'] = df['lag'].map(lambda x: list(np.array(x).cumsum()))
df['cfr_r'] = df['lag_r'].map(lambda x: list(np.array(x).cumsum()))
df.drop(columns=['lag', 'lag_r'], inplace=True)
df = flatten(df).join(data_start, rsuffix='_begin', on='Country/Region')
df['cfr'] = df['cfr'] + df['cfr_begin']
df['cfr_r'] = df['cfr_r'] + df['cfr_begin']
df.drop(columns=['cfr_begin'], inplace=True)
return df.copy()
predict = train.groupby('Country/Region')['lag'].mean()
predict = test[['Date', 'Country/Region', 'lag']].join(predict, on='Country/Region', rsuffix='_r')
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Константный лаг')
ex.add(predict, 'const lag')
import statsmodels.formula.api as smf
from tqdm import tqdm_notebook as tqdm
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = smf.ols('lag ~ n_day1', data=train[train['Country/Region'] == country]).fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.predict(predict[index])
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Линейный лаг')
ex.add(predict, 'linear')
Попробуйте так же применить Double exponential smoothing или Triple exponential smoothing.
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = SimpleExpSmoothing(train[train['Country/Region'] == country].reset_index()['lag']).fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.forecast(sum(index)).values
predict = predict_to_normal(predict, train)
plot_predict(predict, 'SimpleExpSmoothing')
ex.add(predict, 'SimpleExpSmoothing')
from IPython.display import clear_output
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = ExponentialSmoothing(train[train['Country/Region'] == country].reset_index()['lag'], trend='additive').fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.forecast(sum(index)).values
clear_output()
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Double ExpSmoothing')
ex.add(predict, 'Double ExpSmoothing')
Поэксперементируйте с гиперпараметрами модели.
from statsmodels.tsa.arima_model import ARIMA
for param in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (0, 0, param))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA(0, 0, {param}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA(0, 0, {param})')
ex.add(predict, f'ARIMA(0, 0, {param})');
ex.show
Поэксперементируйте с гиперпараметрами модели.
for param in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (param, 0, 0))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA({param}, 0, 0) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA({param}, 0, 0)')
ex.add(predict, f'ARIMA({param}, 0, 0)');
ex.show
Поэксперементируйте с гиперпараметрами модели.
for _AR in range(1, 4):
for _MA in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (_AR, 0, _MA))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA({_AR}, 0, {_MA}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA({_AR}, 0, {_MA})')
ex.add(predict, f'ARIMA({_AR}, 0, {_MA})');
ex.show
_AR = 1
_MA = 1
for _I in range(1, 4):
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = ARIMA(_tr, (_AR, _I, _MA))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index))[0]
except:
bad += 1
if bad:
print(f'ARIMA({_AR}, {_I}, {_MA}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'ARIMA({_AR}, {_I}, {_MA})')
ex.add(predict, f'ARIMA({_AR}, {_I}, {_MA})');
ex.show
from statsmodels.tsa.statespace.sarimax import SARIMAX
_AR = 1
_I = 1
_MA = 1
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
bad = 0
for country in tqdm(all_countries):
try:
_tr = train[train['Country/Region'] == country].reset_index()['lag']
model = SARIMAX(_tr, order=(_AR, _I, _MA))
result = model.fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = result.forecast(steps=sum(index)).values
except:
bad += 1
if bad:
print(f'SARIMAX({_AR}, {_I}, {_MA}) have {bad} errors')
predict = predict_to_normal(predict, train)
plot_predict(predict, f'SARIMAX({_AR}, {_I}, {_MA})')
ex.add(predict, f'SARIMAX({_AR}, {_I}, {_MA})')
Попробуйте использовать библиотеку Prophet для предсказания временного ряда. Документация.
from fbprophet import Prophet
predict = test[['Date', 'n_day1', 'Country/Region', 'lag']].copy()
predict['lag_r'] = np.nan
for country in tqdm(all_countries):
_tr = train[train['Country/Region'] == country].reset_index()[['Date', 'lag']]
_tr.columns = ['ds', 'y']
model = Prophet(weekly_seasonality=True)
result = model.fit(_tr)
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model.predict(model.make_future_dataframe(periods=sum(index)))['yhat'][-sum(index):].values
clear_output()
predict = predict_to_normal(predict, train)
plot_predict(predict, f'Phrofet')
ex.add(predict, f'Phrofet')
Какие признаки оказались наиболее значимыми?
predict = test.copy()
predict['lag_r'] = None
for country in tqdm(all_countries):
model_linear = smf.ols('lag ~ n_day1 + n_day100 + n_day1000 + day_of_week', data=train[train['Country/Region'] == country]).fit()
index = predict['Country/Region'] == country
predict.loc[index, 'lag_r'] = model_linear.predict(predict[index])
if country == 'Russia':
__res = model_linear.summary()
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Линейная модель с регрессией')
ex.add(predict, 'Линейная модель с регрессией')
__res
Наиболее значимые -- число дней после 100 и 1000 заражения, то есть когда ситуация по динамике стабилизировалась.
2.11 (Бонус) Используйте любую другую известную вам модель для предсказания.
from catboost import CatBoostRegressor
def df2normal(df, delay=7, col='lag', for_date=None):
df = df.sort_values('Date').groupby('Country/Region').agg(lambda x: list(x))
delay = datetime.timedelta(delay)
result = []
for country, __df in tqdm(df.iterrows(), total=len(df)):
for l_date, l_nday1, l_nday100, l_nday1000, l_lag in zip(__df['Date'], __df['n_day1'], __df['n_day100'], __df['n_day1000'], __df[col]):
if for_date and for_date != l_date:
continue
row = {
'date': (l_date - train['Date'].min()).days,
'dayofweek': l_date.dayofweek,
'n_day1': l_nday1,
'n_day100': l_nday100,
'n_day1000': l_nday1000,
'lag': l_lag,
'country': country,
}
for r_date, r_nday1, r_nday100, r_nday1000, r_lag in zip(__df['Date'], __df['n_day1'], __df['n_day100'], __df['n_day1000'], __df[col]):
if l_date - delay <= r_date < l_date:
delta = (l_date - r_date).days
row[f'date_{delta}'] = (r_date - train['Date'].min()).days
row[f'dayofweek_{delta}'] = r_date.dayofweek
row[f'n_day1_{delta}'] = r_nday1
row[f'n_day100_{delta}'] = r_nday100
row[f'n_day1000_{delta}'] = r_nday1000
row[f'lag_{delta}'] = r_lag
result.append(row)
result = pd.DataFrame(result)
return result
df3 = df2normal(__df, col='lag')
x_train = df3[df3['date'] < (split_date - train['Date'].min()).days].copy()
x_test = df3[df3['date'] >= (split_date - train['Date'].min()).days].copy()
y_train = x_train['lag']
y_test = x_test['lag']
x_train.drop(columns=['lag'], inplace=True)
x_test.drop(columns=['lag'], inplace=True)
clf_params = {
'depth': 6,
'num_trees': 10000,
'random_strength': 3,
'bagging_temperature': 1.02,
'cat_features': ['country'],
'eval_metric': 'RMSE',
'random_seed': 42,
'use_best_model': True,
'od_type': 'Iter',
'od_wait': 200,
'task_type': 'CPU',
}
clf = CatBoostRegressor(**clf_params)
clf.fit(x_train, y_train, eval_set=(x_test, y_test), plot=True, verbose=False)
clf.get_feature_importance(prettified=True)
df4 = __df.copy()
df4['lag'] = np.nan
for date in x_test['date'].unique():
date = train['Date'].min() + datetime.timedelta(int(date))
df4_ = df2normal(df4, col='lag', for_date=date)
df4_.drop(columns='lag', inplace=True)
df4.loc[df4['Date'] == date, 'lag'] = clf.predict(df4_)
predict = test.copy()
predict['lag_r'] = df4[~df4['lag'].isna()]['lag']
predict = predict_to_normal(predict, train)
plot_predict(predict, 'Catboost')
ex.add(predict, 'Catboost')
ex.show.to_pickle('exp_results.pkl')
models = {k: 0 for k in ex.show.index}
top_models = {k: 0 for k in ex.show.index}
for col in ex.show.columns:
for i, model in enumerate(ex.show.sort_values(col).index):
if i == 0:
top_models[model] += 1
models[model] += i
Топ моделей по количеству первых мест:
pd.Series(top_models).sort_values(ascending=False)
Рейтинг моделей по итоговой позиции:
pd.Series(models).sort_values(ascending=True)
Как можно заметить, catboost победил все модели по качеству предсказания (выборка делилась по времени, данные для предсказания каждый раз обновлялись на основе предсказаний).
Тем не менее, некоторые алгоритмы близко подходят по качеству. Возможно, catboost выиграл за счет обобщающей способности по всем странам, возможно ряд не совсем классический для задач временных рядов.
Кроме того можно провести дополнительные эксперименты по предсказанию для catboost на cfr, а не лаги; попробовать применить LSTM и другие RNN.
Тем не менее, для каждой задачи возможно лучше будет использовать свой алгоритм, который лучше всего описывает данный ряд.